
#include "functions.h"
#include "mrandom.hpp"
#include <cmath>
#include <fstream>
#include <iostream>
#include <string>
using namespace std;
using namespace Eigen;

constexpr int u{0}, v{1};

void ReadMeshASSCII(const std::string &meshFname, VarArray<Point, 2, 0> &grid) {
  cout << "start reading mesh file\n";

  ifstream meshF{meshFname};
  if (!meshF) {
    cerr << "open mesh file err\n" << endl;
    exit(EXIT_FAILURE);
  }

  int ni, nj;
  meshF >> ni >> nj;

  // 跟 blazik 同步 后面有需求再改
  grid.ReConstruct(ni, nj);

  for (int j = 0; j < nj; ++j)
    for (int i = 0; i < ni; ++i) {
      meshF >> grid(i, j).x;
    }

  for (int j = 0; j < nj; ++j)
    for (int i = 0; i < ni; ++i) {
      meshF >> grid(i, j).y;
    }

  meshF.close();
  cout << "read mesh file success" << endl;
}
/**
void InitializeFlow(VarArray<Variable, 2, 0> &flow) {
  const int np_in_ydir = flow.Dim(1);
  flow(0, 0).u = 0.0;
  flow(0, 0).v = 0.0;
  // 给第一列赋初值
  for (int j = 1; j < np_in_ydir; ++j) {
    flow(0, j).u = 1.0;
    flow(0, j).v = 0.0;
  }
  const int np_in_xdir{flow.Dim(0)};
  // 给上下边界赋值
  for (int i = 1; i < np_in_xdir; ++i) {
    flow(i, 0).u = 0.0;
    flow(i, 0).v = 0.0;

    flow(i, np_in_ydir - 1).u = 1.0;
    flow(i, np_in_ydir - 1).v = 0.0;
  }
}
**/
void InitializeFlow(VarArray<Vector2d, 2, 0> &flow) {
  const int np_in_ydir = flow.Dim(1);
  Vector2d &var00 = flow(0, 0);
  var00(0) = 0.0;
  var00(1) = 0.0;
  // 给第一列赋初值
  for (int j = 1; j < np_in_ydir; ++j) {
    Vector2d &var0j{flow(0, j)};
    var0j(0) = 1.0;
    var0j(1) = 0.0;
  }
  const int np_in_xdir{flow.Dim(0)};
  // 给上下边界赋值
  for (int i = 1; i < np_in_xdir; ++i) {
    Vector2d &vari0{flow(i, 0)};
    vari0(0) = 0.0;
    vari0(1) = 0.0;

    Vector2d &varIlast{flow(i, np_in_ydir - 1)};
    varIlast(0) = 1.0;
    varIlast(1) = 0.0;
  }
}

/**
// 计算之前未知时间层需要使用已知时间层初始化
void CalcJacobianMatrix(VarArray<Variable, 2, 0> &flow, int column,
                        const VarArray<Point, 2, 0> &grid,
                        VarArray<Matrix2d, 2, 0> &jacob_matrix,
                        const double re_num) {
  int npoint_in_ydir{grid.Dim(1)};
  const double C{1 / re_num};

  // 先统一处理 最后考虑首尾
  for (int j = 1; j < npoint_in_ydir - 1; ++j) {

    Matrix2d &matJm1{jacob_matrix(0, j)};
    Matrix2d &matJ{jacob_matrix(1, j)};
    Matrix2d &matJp1{jacob_matrix(2, j)};

    const double rJacob = 2.0 / (grid(column, j + 1).y - grid(column, j - 1).y);
    const double B =
        grid(column, j + 1).y - 2 * grid(column, j).y + grid(column, j - 1).y;
    const double rDeltx = 1.0 / (grid(column, j).x - grid(column - 1, j).x);
    // 第一次迭代第i层时间步上的值由上一层近视
    // if(iter_num != 1)

    matJm1(0, 0) = 0.0;
    matJm1(1, 0) = -0.5 * flow(column, j).v * rJacob - C * rJacob * rJacob -
                   0.5 * B * rJacob * rJacob * C * rJacob;
    matJm1(0, 1) = 0.5 * rJacob;
    matJm1(0, 0) = 0.0;

    matJ(0, 0) = 1.0 * rDeltx;
    matJ(1, 0) = (2.0 * flow(column, j).u - flow(column - 1, j).u) * rDeltx +
                 2.0 * C * rJacob * rJacob;
    matJ(0, 1) = 0.0;
    matJ(0, 0) = 0.5 * (flow(column, j + 1).u - flow(column, j - 1).u) * rJacob;

    matJ(0, 0) = 0.0;
    matJ(1, 0) = 0.5 * flow(column, j).v * rJacob - C * rJacob * rJacob +
                 0.5 * B * rJacob * rJacob * C * rJacob;
    matJ(0, 1) = 0.5 * rJacob;
    matJ(0, 0) = 0.0;
  }
  // 首尾处理
  jacob_matrix(0, 1) << 0.0, 0.0, 0.0, 0.0;
  jacob_matrix(2, npoint_in_ydir - 2) << 0.0, 0.0, 0.0, 0.0;
}
**/

// /**
// 计算之前未知时间层需要使用已知时间层初始化
void CalcJacobianMatrix(const VarArray<Vector2d, 2, 0> &flow, int column,
                        const VarArray<Point, 2, 0> &grid,
                        VarArray<Matrix2d, 2, 0> &jacob_matrix,
                        const double re_num) {
  int npoint_in_ydir{grid.Dim(1)};
  const double C{1 / re_num};

  // 先统一处理 最后考虑首尾
  for (int j = 1; j < npoint_in_ydir - 1; ++j) {

    Matrix2d &matJm1{jacob_matrix(0, j)};
    Matrix2d &matJ{jacob_matrix(1, j)};
    Matrix2d &matJp1{jacob_matrix(2, j)};

    const Vector2d &varIJ{flow(column, j)}, varIm1J{flow(column - 1, j)},
        varIJm1{flow(column, j - 1)}, varIJp1{flow(column, j + 1)};

    const double rJacob = 2.0 / (grid(column, j + 1).y - grid(column, j - 1).y);
    const double B =
        grid(column, j + 1).y - 2 * grid(column, j).y + grid(column, j - 1).y;
    const double rDeltx = 1.0 / (grid(column, j).x - grid(column - 1, j).x);

    // cout << rJacob << "\t " << B << "\t " << rDeltx << endl;
    // 第一次迭代第i层时间步上的值由上一层近视
    // if(iter_num != 1)
    matJm1(0, 0) = 0.0;
    matJm1(1, 0) = -0.5 * varIJ(1) * rJacob - C * rJacob * rJacob -
                   0.5 * B * rJacob * rJacob * C * rJacob;
    matJm1(0, 1) = 0.5 * rJacob;
    matJm1(1, 1) = 0.0;

    matJ(0, 0) = 1.0 * rDeltx;
    matJ(1, 0) =
        (2.0 * varIJ(0) - varIm1J(0)) * rDeltx + 2.0 * C * rJacob * rJacob;
    matJ(0, 1) = 0.0;
    matJ(1, 1) = 0.5 * (varIJp1(0) - varIJm1(0)) * rJacob;
    // cout << ma
    matJp1(0, 0) = 0.0;
    matJp1(1, 0) = 0.5 * varIJ(1) * rJacob - C * rJacob * rJacob +
                   0.5 * B * rJacob * rJacob * C * rJacob;
    matJp1(0, 1) = 0.5 * rJacob;
    matJp1(1, 1) = 0.0;
  }
  // 首尾处理
  jacob_matrix(0, 1) << 0.0, 0.0, 0.0, 0.0;
  jacob_matrix(2, npoint_in_ydir - 2) << 0.0, 0.0, 0.0, 0.0;
}
// **/
void CalcFXk(VarArray<Vector2d, 2, 0> &flow, int column,
             const VarArray<Point, 2, 0> &grid, const double re_num,
             vector<Vector2d> &fx_k)

{
  int npoint_in_ydir{grid.Dim(1)};
  const double C{1 / re_num};

  for (int j = 1; j < npoint_in_ydir - 1; ++j) {
    const double rJacob = 2.0 / (grid(column, j + 1).y - grid(column, j - 1).y);
    const double B =
        grid(column, j + 1).y - 2 * grid(column, j).y + grid(column, j - 1).y;
    const double rDeltx = 1.0 / (grid(column, j).x - grid(column - 1, j).x);

    const Vector2d &var_i_j{flow(column, j)}, var_im1_j{flow(column - 1, j)},
        var_i_jm1{flow(column, j - 1)}, var_i_jp1{flow(column, j + 1)};

    Vector2d &fx_j{fx_k[j]};
    fx_j(0) = (var_i_j(u) - var_im1_j(u)) * rDeltx +
              0.5 * (var_i_jp1(v) - var_i_jm1(v)) * rJacob;

    fx_j(1) = var_i_j(u) * (var_i_j(u) - var_im1_j(u)) * rDeltx +
              0.5 * var_i_j(v) * (var_i_jp1(u) - var_i_jm1(u)) * rJacob -
              C * rJacob * rJacob *
                  ((var_i_jp1(u) - 2 * var_i_j(u) + var_i_jm1(u)) -
                   0.5 * (var_i_jp1(u) - var_i_jm1(u)) * rJacob * B);
    // fx_k(j)(0) = 0.0;
    // cout << fx_j << endl;
  }
}

double CalcDifference(const vector<Vector2d> &vec1,
                      const vector<Vector2d> &vec2);
/**
// 高斯晒德尔迭代 计算残差
void CalcResidual(const VarArray<Matrix2d, 2, 0> &jacob_matrix,
                  const vector<Vector2d> &fx_k, vector<Vector2d> &residual,
                  const VarArray<Point, 2, 0> &grid) {
  Vector2d zero;
  zero << 0.0, 0.0;
  for (auto &vec : residual)
    vec = zero;
  // cout << residual[j] << endl;

  int npoint_in_ydir = residual.size();
  vector<Vector2d> residual_old = residual;

  int count = 1;
  Vector2d mid;
  double res1, resn, tmp;
  do {
    for (int j = 1; j < npoint_in_ydir - 1; ++j) {

      // if (fabs(jacob_matrix(1, j).determinant()) < 1.0e-16) {
      //   cout << j << endl;
      //   cout << (jacob_matrix(1, j)) << endl;
      //   cerr << "cannot get inverse \n" << endl;
      //   exit(EXIT_FAILURE);
      // }

      // residual[j] = jacob_matrix(1, j).inverse() *
      //               (fx_k[j] - jacob_matrix(0, j) * residual[j - 1] -
      //                jacob_matrix(2, j) * residual[j + 1]);
      Vector2d &resJ{residual[j]};
      mid = (fx_k[j] - jacob_matrix(0, j) * residual[j - 1] -
             jacob_matrix(2, j) * residual[j + 1]);

      // const double deltx = (grid(1, j).x - grid(0, j).x);
      const double deltx = 1.0e-7;
      const Matrix2d &mat{jacob_matrix(1, j)};
      resJ[0] = mid(0) * deltx;
      resJ[1] = (mid(1) - resJ[0] * mat(1, 0)) / mat(1, 1);
      // cout << jacob_matrix(1, j).inverse() << endl;
      // cout << residual[j] << endl;
    }
    //
    resn = CalcDifference(residual, residual_old);
    residual_old = residual;

    if (count == 1)
      res1 = resn;
    ++count;
    tmp = log10(resn / res1);
    cout << "tmp = " << tmp << endl;
    cout << "resn = " << resn << endl;
  } while (tmp > -8.0 && resn > 1.0e-8 && count < 80);

  cout << "count = " << count << endl;
}
**/
vector<double> CalcResidual2(const VarArray<Matrix2d, 2, 0> &jacob_matrix,
                             const vector<Vector2d> &fx_k) {
  int npoint_in_ydir = jacob_matrix.Dim(1);

  VarArray<double, 2, 0> coeff(3, npoint_in_ydir - 2);
  vector<double> diff_u(npoint_in_ydir - 2);

  for (int i = 0; i < npoint_in_ydir - 2; ++i) {

    const Matrix2d &mat1 = jacob_matrix(0, i + 1);
    coeff(0, i) = mat1(1, 0);
    const Matrix2d &mat2 = jacob_matrix(1, i + 1);
    coeff(1, i) = mat2(1, 0);
    const Matrix2d &mat3 = jacob_matrix(2, i + 1);
    coeff(2, i) = mat3(1, 0);
    const Vector2d &fx = fx_k[i + 1];
    diff_u[i] = fx(1);
  }
  ThomasAlgorithm(coeff, diff_u.data(), 0, diff_u.size());
  return diff_u;
}

double CalcDifference(const vector<Vector2d> &vec1,
                      const vector<Vector2d> &vec2) {
  double difference{0.0};
  int size = vec1.size();
  for (int i = 1; i < size - 1; ++i) {
    const Vector2d &v1{vec1[i]}, v2{vec2[i]};
    difference +=
        (v1(u) - v2(u)) * (v1(u) - v2(u)) + (v1(v) - v2(v)) * (v1(v) - v2(v));
  }
  difference = sqrt(difference);
  return difference;
}

void UpdateFlow(VarArray<Vector2d, 2, 0> &flow,
                const vector<Vector2d> &residual, int column) {

  int npoint_in_ydir = residual.size();

  for (int j = 1; j < npoint_in_ydir - 1; ++j) {
    flow(column, j) -= residual[j];
  }
}

double CalcNormal(const vector<Vector2d> &residual) {
  int size = residual.size();
  double tmp = 0.0;
  for (int i = 1; i < size - 1; ++i) {
    const Vector2d &vec{residual[i]};
    tmp += vec(u) * vec(u) + vec(v) * vec(v);
  }
  return sqrt(tmp);
}

// 有些问题
void InitializeNextLine(VarArray<Vector2d, 2, 0> &flow, int column) {
  int npoint_in_ydir{flow.Dim(1)};
  int npoint_in_xdir{flow.Dim(0)};
  double delta = 1.0 / (npoint_in_ydir - 1);
  // Rand_double randnum(0.0, 1.0e-6);
  // Vector2d tmp;
  for (int j = 1; j < npoint_in_ydir - 1; ++j) {
    // tmp(0) = randnum();
    // tmp(1) = randnum();
    // flow(column, j) = flow(column - 1, j) + tmp;
    Vector2d &var_i_j = flow(column, j);
    var_i_j(u) = j * delta;
    var_i_j(v) = j * delta * 0.1;
  }
}

template <typename T> void PrintForCheck(const VarArray<T, 2, 0> &vv) {
  int ni = vv.Dim(0);
  int nj = vv.Dim(1);

  for (int j = 0; j < nj; ++j) {
    for (int i = 0; i < ni; ++i) {
      cout << (vv(i, j)) << "\t";
    }
    cout << endl << endl;
  }
}